Calculation of nuclear matrix elements in neutrinoless double electron capture 
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We compute nuclear matrix elements for neutrinoless double electron capture on 152 Gd, 164 Er 
and 180 W nuclei. Recent precise mass measurements for these nuclei have shown a large resonance 
enhancement factor that makes them the most promising candidates for observing this decay mode. 
We use an advanced energy density functional method which includes beyond mean-field effects such 
as symmetry restoration and shape mixing. Our calculations reproduce experimental charge radii 
and B(E2) values predicting a large deformation for all these nuclei. This fact reduces significantly 
the values of the NMEs leading to half-lives larger than 10 29 years for the three candidates. 
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I. INTRODUCTION 

The nature of the neutrino as a Majorana or Dirac 
particle is still an open question. Experimental evidence 
of lepton-number violation weak processes such as neu- 
trinoless double beta decay (0^/3/3) or double electron 
capture (O^ECEC) would determine unambiguously the 
Majorana character of this elementary particle. In par- 
ticular, in O^ECEC a Majorana neutrino is exchanged in 
the capture of two electrons (e~) from the inner shells of 
an atom with mass number A and number of protons Z: 



(A,Z)^(A,Z-2)* 



(1) 



The daughter atom (A, Z — 2)** is unstable and decays 
to its ground state (A, Z — 2) by X-rays and/or Auger 
electrons emission from the electronic excitation and by 
7-ray emission if the nucleus is in an excited state. In 
Fig. 1 we show a schematic view of the process and its 
energetics. The capture rate Xouee can be written as [1]: 
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where A = Q ee — £?**, Q ee is the difference between 
the initial and final atomic masses, E** and T^ are the 
energy and the width of the excited state of the daugh- 
ter atom with two electronic holes in states a, b and an 
excited nucleus, and V a b is the transition amplitude be- 
tween the initial and final states [2] : 
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(3) 



with ra/3/3, Gp: Oc, Qa = 1-25, R being the effective Ma- 
jorana neutrino mass, the Fermi constant, the Cabbibo 
angle, the axial- vector coupling constant and the nuclear 
radius respectively; (F a b) is the lepton part of the ma- 
trix element that takes into account the atomic electron 
wave functions of the states a, b and M 0l/ is the nuclear 
matrix element (NME) [2] . Equation 2 shows a resonant 
behavior with an increase of the capture rate whenever 
A ~ T is fulfilled, otherwise the process is very much sup- 
pressed. This is a very restricted condition to search for 
isotopes with a potential O^ECEC because typical values 



for the widths of the atomic states are of the order of 
electron volts (eV) while Q ee is found by subtracting two 
quantities which are of the order of MeV. Recent experi- 
ments using Penning traps have reported extremely pre- 
cise measurements of Q ee values for the most promising 
O^ECEC candidates where the initial and final nuclei are 
in their ground states, namely, 152 Gd [3], 164 Er [4] and 
180 W [5]. These experiments have shown rather large res- 
onance enhancement factors that could lead to O^ECEC 
half-lives between 10 26-30 years/ {rnpp[eV]) 2 [5]. How- 
ever, this half-lives are based on NMEs computed us- 
ing spherical quasiparticle random phase approximation 
(QRPA) for mother and daughter nuclei [2]. Recently 
Fang et al. [6] have extended these calculations to de- 
formed shapes and found a substantial suppression of the 
NMEs that leads to longer half-lives. However, in these 
calculations the deformation is a free parameter that is 
not consistently determined within the model. 

In this work we use an energy density functional (EDF) 
method including beyond-mean field correlations to eval- 
uate the O^ECEC NMEs for these isotopes. This method 
was originally developed to study in a self-consistent 
manner the role of deformation, shape mixing and pair- 
ing correlations in 0uf3f3 NMEs [7, 8]. Here, a density- 
dependent effective two body interaction between the 
constituent nucleons is assumed as the starting point. 
Then, single particle energies, residual interactions be- 
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FIG. 1. Energetics of the Oz^ECEC process. See text for 
details. 



yond the mean field, pairing and deformation content 
of the nuclei, transitions and decays rates, etc., are de- 
duced self-consistently from solving the associated quan- 
tum many-body problem. In addition, mixing of differ- 
ent shapes is also allowed within this formalism. This is a 
different perspective that the one used in QRPA calcula- 
tions [6] , where the several free parameters of the model 
are adjusted specifically for each capture. In particular, 
the deformation of initial and final nuclei are fitted to 
reproduce either the spectroscopic quadrupole moment 
of the 2+ state or from the reduced transition probabili- 
ties B(E2). The effects of shape mixing, particle number 
and angular momentum restoration are not considered. 
In addition, pairing gaps are also found by fitting at the 
BCS level experimental even-odd mass differences. Fi- 
nally, single particle energies and residual interactions do 
not come from the same underlying interaction, and an 
additional renormalization is required adding two more 
adjustable parameters to the model. 



II. THEORETICAL FRAMEWORK 

We now summarize the formalism that we use to com- 
pute the NMEs considering that the initial and final nu- 
clei are in their ground states. First of all, M 0u can be 
expressed as the sum of the so-called Fermi (F), Gamow- 
Teller (GT) and tensor (T) parts [9]: 



M 0u = 
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with gy = 1. In this work we neglect the tensor term that 
it is estimated to be small [10, 11]. The other terms in Eq. 
4 are written as the overlaps of two-body operators be- 
tween initial and final states, M^f jGT = (0^|M^ GT |0+), 
with: 



r(i)r(2) 
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Mp = V F ^>, M^ T = V GT (&^-&^)i^i^ (5) 

Here i+ is the isospin ladder operator that changes pro- 
tons into neutrons and a are the Pauli matrices acting on 
the spin part. The spatial part of the wave functions is 
affected by the so-called neutrino potentials Vf/gt that 
include nucleon finite size corrections [12], radial short- 
range correlations within the unitary correlation operator 
Method (UCOM) [13, 14] and higher order currents [12]. 
Detailed expressions for these operators can be found in 
Refs. [10, 12]. 

Initial (i) and final (/) states are found by using config- 
uration mixing techniques based on the generator coor- 
dinate method (GCM) with particle number and angular 
momentum restoration of Hartree-Fock-Bogoliubov type 
(HFB) mean-field wave functions [15-18]: 
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where the axial quadrupole deformation fa is used as the 
generating coordinate and gLAfa) are the coefficients 



of the linear combination which are obtained by solv- 
ing the corresponding Hill- Wheeler- Griffin (HWG) equa- 
tions [18]. The set of wave functions Y^hAfa)) are de- 
fined as: 



|M/f // (/? 2 ))=P / P iY «//P^//|$ l// (/? 2 )) 



(7) 



being P 7 , P N and P z the projection operators on to 
angular momentum / and number of neutrons and pro- 
tons, respectively; \<&i/f(fa)) are HFB type mean- field 
wave functions that are found by using the variation af- 
ter particle number projection method [19]. We have 
used two parametrizations of the Gogny force as the un- 
derlying density-dependent interaction, namely, the most 
widely used in nuclear structure calculations D1S [20] 
and the most recent one DIM [21] that is designed to 
give the best fit to experimental masses in the whole 
nuclear chart. Experimental observables such as nu- 
clear masses, radii, spectra, transition rates, etc. can 
be directly compared to the expectation values calcu- 
lated from the GCM wave functions given in Eq. 6. In 
addition, the overlaps between the projected states de- 
fined in Eq. 7 of the corresponding operators give the 
dependence on the quadrupole deformation of such ob- 
servables. Finally, probability distributions for the GCM 
wave functions -Eq. 6- to have a given deformation, the 
so-called collective wave functions FhAfa), provide ad- 
ditional information about the intrinsic structure of these 
states and identify the relevant regions within the space 
parametrized by the deformation fa (see [7, 8] and refer- 
ences therein for detailed expressions). In this work only 
HFB axial symmetric configurations are considered and 
these intrinsic wave functions are expanded in a single 
particle basis composed by eleven major spherical har- 
monic oscillator shells. 



III. RESULTS 



We now present the results obtained for the O^ECEC 
processes 152 Gd -+ 152 Sm, 164 Er -> 164 Dy and 180 W -+ 
180 Hf. First, to check the reliability of the method ap- 
plied to these nuclei, we compare in Table I theoretical 
and experimental values for different observables such as: 
binding energy per particle and nuclear charge radius - 
ground state properties- or excitation energy of the first 
2+ state and its electric quadrupole reduced transition 
probability B(E2;2^ — » 0^) -spectroscopic properties. 
On the one hand, we see that the results provided by 
D1S and DIM parametrizations are almost identical, al- 
though D1S gives a bit higher binding energies, radii and 
B(E2) values. On the other hand, theoretical predic- 
tions give a bit smaller binding energies and radii, larger 
B(E2) values and, except for A = 152, larger 2 + excita- 
tion energies than the experimental data. Nevertheless, 
the agreement with the experiments is very good, taking 
into account that the interactions are globally fitted and 
effective charges are not needed within this framework. 



Furthermore, the relative systematics for the pair of nu- 
clei involved in the electron capture is reproduced, i.e., 
whenever a quantity is experimentally larger or smaller in 
one nucleus than in the other, the same trend is observed 
in the theoretical results. Other degrees of freedom not 
included in these calculations like, for example, triaxial- 
ity and time-reversal symmetry breaking could help to 
improve the agreement with the experimental data, es- 
pecially for the excitation energies and transition proba- 
bilities. Work in this direction is in progress. 

All the isotopes studied in this work are open shell 
nuclei, both in protons and neutrons. This means that 
collective effects are expected to appear. In fact, most 
of these isotopes show experimentally rotational ground 
state bands [22] pointing out that quadrupole deforma- 
tions play a significant role. The GCM method used here 
is particularly suitable to include this degree of freedom 
in a self-consistent manner. In Figs. 2(a)- 3(a)- 4(a) we 
plot the collective wave functions -distribution of proba- 
bility of having the state a given deformation- as a func- 
tion of the quadrupole deformation fa for the nuclei cal- 
culated in this work using the D1S parametrization (sim- 
ilar features are found for DIM). We observe that all the 
isotopes have a large prolate deformation. Smaller defor- 
mation and narrower distribution is obtained for 152 Gd 
than for 152 Sm. The maxima of the collective wave func- 
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D1S 


8.212 


5.027 


0.227 


113 


152 Gd 


DIM 


8.181 


4.991 


0.228 


104 




Expt 


8.233 


5.082 


0.344 


73(7) 




D1S 


8.219 


5.061 


0.118 


204 


152 Sm 


DIM 


8.195 


5.024 


0.117 


186 




Expt 


8.244 


5.084 


0.122 


144(3) 




D1S 


8.118 


5.190 


0.125 


226 


164 Er 


DIM 


8.088 


5.154 


0.125 


215 




Expt 


8.149 


5.238 


0.091 


218(7) 




D1S 


8.122 


5.181 


0.105 


231 


1M Dy 


DIM 


8.103 


5.149 


0.103 


225 




Expt 


8.158 


5.221 


0.073 


209(3) 




D1S 


7.982 


5.321 


0.137 


172 


180 W 


DIM 


7.955 


5.289 


0.136 


172 




Expt 


8.025 


5.349 


0.093 


- 




D1S 


7.987 


5.312 


0.133 


175 


180 Hf 


DIM 


7.968 


5.276 


0.130 


167 




Expt 


8.034 


5.342 


0.104 


155(5) 



tions are found at fa ~ 0.25 and fa ~ 0.35 respectively. 
For the other two cases mother and daughter collective 
wave functions are very similar, having their maxima at 
fa ~ 0.35 for A = 164 and fa ~ 0.30 for A = 180. There- 
fore, the quadrupole degree of freedom has to be taken 
into account in the calculation of the O^ECEC rates. 

The dependence of the NMEs on the deformation can 
be obtained by calculating the overlaps of the operators 
defined in Eq. 5 between the particle number and an- 
gular momentum projected states given in the expres- 
sion 7 [7]. We represent the intensity of the Gamow- 
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TABLE I. Theoretical and experimental data for ground 
state and spectroscopic properties of the nuclei involved in 
Oz^ECEC. First and second rows of each cell correspond 
to theoretical values calculated with Gogny D1S and DIM 
parametrizations respectively. Experimental data are taken 
from Refs. [22-24] 



FIG. 2. (color online) (a) Collective wave functions ; (b) in- 
tensity of the Gamow- Teller part of the Oz/ECEC NME and 
(c) pairing energies as a function of the quadrupole defor- 
mation fa for 152 Gd and 152 Sm isotopes calculated with the 
Gogny D1S parametrization. 



P 2 



-0.4 -0.2 



(a) GOGNYD1S 

— 164 Er (0, + ) 
-- 164 Dy(0;) 




-0.4 



-0.2 



P 2 
0.2 



0.4 



0.6 




(a) GOGNYD1S /\ 
— 180 W (0^) 

-— 180 Hf(0 f + ) 





FIG. 3. (color online) Same as Fig. 2 but for 164 Er and 164 Dy 
isotopes. 



FIG. 4. (color online) Same as Fig. 2 but for 180 W and 180 Hf 
isotopes. 



Teller part of the NME as a function of the deformation of 
the mother and daughter nuclei in Figs. 2(b)- 3(b)- 4(b). 
The same qualitative behavior is found for the Fermi 
part and also for the DIM parametrization and they are 
not shown. The three cases display similar characteris- 
tics. Highest intensities are located around the spher- 
ical shape and along the diagonal ft™ 1 = /3§ n between 
fa € [-0.3,0.25] for A = 152,164 and /3 2 G [-0.2,0.2] 
for A = 180. The rest of configurations are very much 
suppressed. Therefore, O^ECEC is favored whenever the 
mother and daughter are nearly spherical in these cases. 
We also observe that the distribution of the NMEs as 
a function of the deformation is correlated to the pair- 
ing energy (E pp ) of the nuclei involved in the capture 



process (see Figs. 2(c)- 3(c)- 4(c)). Hence, whenever a 
maximum in the paring energy is present, also a maxi- 
mum in the NME is obtained. The same correlation has 
been accounted for in the 0vj3f3 NMEs [7, 8] and using 
a seniority scheme [14, 25]. However, although spherical 
shapes are favored by the Oz/ECEC operator, we see in 
Figs. 2(b)- 3(b)- 4(b) that the actual many-body wave 
functions explore only regions that corresponds to quite 
small intensities- see the shaded area. Therefore, the 
NMEs are strongly suppressed whenever the deforma- 
tion of the nuclei is taken into account as it was already 
observed in the calculation of 0vf3f3 NMEs of 150 Nd [7]. 
A similar qualitative result is obtained in QRPA calcu- 
lations [6] although the suppression with respect to the 



Nucleus 


Par am. 


EDF 


T^ n (y) 


M 0u 
QRPA [6] 


152 Gd 


DIS 
DIM 


1.07 
0.89 


4.2 x 10 29 
6.2 x 10 29 


3.23-2.67 


164 Er 


DIS 
DIM 


0.64 
0.50 


1.3 x 10 34 
2.1 x 10 34 


2.64-2.27 


180 W 


DIS 
DIM 


0.58 
0.38 


1.6 x 10 32 
3.8 x 10 32 


2.05-1.79 



TABLE II. NMEs and estimated half- lives for most probable 
Ovee captures. 



spherical result is much larger in this work. 

To shed more light on the structure of the NMEs and 
pairing energies shown above, we plot in Fig. 5 the single 
particle energies (s.p.e.) as a function of the deforma- 
tion for the A = 164 isotopes. A similar analysis applies 
also to the rest of isotopes studied here. In the spheri- 
cal point (/^2 = 0) these levels have a degeneracy equal 
to (2j + 1). For /? 2 7^ this degeneracy is broken in 
(j + 1/2) doubly degenerated levels. In the following dis- 
cussion we maintain the name of the spherical orbits for 
referring the different levels. It is well known that pair- 
ing energies are larger whenever the level density around 
the Fermi energy is larger. This is precisely what it is ob- 
served in the range of deformations 02 £ [—0.2, 0.2] where 
proton levels coming from the orbits 7^3/2, ^11/2? ^l^ 
are very close in energy. For the neutrons, the fact that 
the spherical orbit uh 9 /2 is not fully occupied, and the 
fast lowering of the ^3/2, ^13/2 levels and a rising of the 
vfj/2 induced by prolate deformation leads to a large 
level density around spherical shapes. Hence, large pair- 
ing correlations are expected there and, correspondingly, 
large NMEs are predicted. On the other hand, for larger 
prolate deformations {fa £ [0.3,0.5]) we observe that the 
density of levels around the Fermi level is smaller for pro- 
tons and neutrons. In particular, gaps are found in the 
proton s.p.e. produced by the emptying of some levels 
coming from the 71^3/2, ^5/2 , ^97/2 an d the filling of in- 
truder nhn/2 orbits. In the neutron part, the reduction 
of the level density is due to the sharp decrease of some 
^13/2 levels, which start to be occupied at fa ~ 0.2, and 
the increase of some orbits from ^9/2, vf*j y 2 . 

The final results for the NMEs and half-lives calculated 
with DIS and DIM parametrizations are summarized in 
the Table II. We observe smaller values of the NMEs 
for the DIM parametrization although the difference be- 
tween the results is not very significant. The largest NME 
corresponds to the decay of 152 Gd that, together with a 
large resonance enhancement, results in the shortest half- 
live among the candidates studied here (T1/2 ~ 10 29 y). 
In the evaluation of the half-lives we assume a neutrino 
Majorana mass mpp = 50 meV and the values of (Fab), A 
and T a 5 are taken from Ref. [6] corresponding to the elec- 
tronic transition that gives the shorter half-life in each 



> 

CD 



Q- 



'/~<S-rb 


\\\ \ /// // 1 


1 V 1 / l VV*. ' 




164 Er ^ 


0% \ n /" / 1 


y (a) neutrons '1 


il64 Dy 


/ \ V /\ 






NN ^^C^o^io J' ' 


\ /£ 






\ / \\ 


\ \ / 


V./p3fc N X\ 


x \ > 


r\ \ / 


// '%\ 




Y ^*~y 




\" v/ ,! 






*3( / 




v\ -• y\ 


/ \ / ' 
/ V V 


C- - - J- - > %" -VN- 


\ \ V \ \ 
\\ ^ \ / 


/'vf 


/ / \ \ N 




/ / A 


/ / \ \\ 


W 1 \ -A/ 


' y/v 


' / *~^ \ \ 


\ \V^ \ /. ^^ 


- f '/ A 


/ / f\ \ V 


\ A v/ \ - 


v> A ' 


(82) \ > 


v A \ //\ A 


W/ / \ ' 


V — ' \ \ 


\ ' \ ^y/ \ a 


X \' 


\ 


v \ / \ , / Y 


\\ \ 1 \t.\ 


1 1 1 1 1 1 1 l\ 


\ v 1 1 \ 1 / 1 1 /K / 




FIG. 5. (color online) Single particle energies (continuous 
and dashed lines represent positive and negative parity levels 
respectively) as a function of the deformation for A — 164 for 
(a) neutrons and (b) protons. The red diamonds and the blue 



bullets are the Fermi levels for ^Erge and 
respectively. 



gDygs isotopes 



case. We also compare the NMEs obtained in this work 
with the most recent deformed QRPA calculations [6]. 
Our values are a factor ~ 3 — 4 smaller than the QRPA 
ones. However, larger NMEs are obtained within our ap- 
proach if we neglect the configuration mixing and assume 
single deformations for initial and final states similar to 
the ones used in QRPA calculations. This shows again 
the relevance of the deformation and configuration mix- 
ing in determining the value of the NMEs. 



IV. CONCLUSSIONS 

In summary, we have used a state-of-the-art energy 
density functional method to compute nuclear matrix 
elements for the most promising O^ECEC candidates. 



Our approach includes particle number and angular mo- 
mentum symmetry restoration and shape mixing along 
the axial quadrupole deformation /^2, using Gogny D1S 
and DIM as the underlying interactions. We observe a 
strong reduction of the NMEs whenever the deformation 
of mother and daughter nuclei is not close to the spherical 
shape. All isotopes studied here have a large prolate de- 
formation, and therefore, we predict very small values of 
the NMEs. Finally, for all cases considered in this work 
we find half-lives longer than 10 29 years despite the large 
resonance enhancement found in very recent experiments 



for these nuclei. 
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